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Abstract. The fragmented population of alpine salamanders from the Dinaric Alps belongs to a distinct evolutionary line- 
age known as Salamandra atra prenjensis. However, the phylogenetic relationships within this lineage are unknown since 
previous studies did not comprehensively include Dinaric population fragments. In this study, we use six microsatellite 
loci and two mtDNA markers to examine the evolutionary relationships and genetic structure in several isolated fragments 
of alpine salamanders that are well spread along the Dinarides. We discovered that during Pleistocene glaciations, S. atra 
persisted in at least two Dinaric refugia: an old one in the Prenj mountain, and a more recent one in Gorski Kotar, whereas 
there is evidence of colonization for the populations of the Cvrsnica and Prokletije mountains. The results revealed that 
mountain Prenj was probably the main diversification center of alpine salamanders in the Dinarides. In addition, to pro- 
vide a state-of-the-art review of the evolutionary history of the species, we also included available sequences of S. atra from 
the entire distribution range; we discuss the results from a conservation perspective. 
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Introduction 


In the current epoch of the sixth mass extinction, charac- 
terized by a widespread and rapid decline of biodiversity 
(RIPPLE et al. 2017), phylogenetic, phylogeographic and 
population genetic studies are becoming increasingly im- 
portant in the context of conservation biology and species 
management (FRANKHAM et al. 2002). 

The current geographical patterns of genetic variation 
in many species are a result of climatic and geomorpho- 
logical changes that happened throughout the Earth’s ev- 
olution (LOMOLINO & HEANEY 2004). Probably the most 
famous example of such events are the climatic oscillations 
in the Pleistocene that caused a series of alternated con- 
tractions and expansions of distributional ranges (HEWITT 
1999, 2000, 2004). Traditionally, European species were 
thought to have survived the Pleistocene glaciations in a 


few large, continuous areas in the extreme south of Europe 
- an idea largely bolstered by paleoclimatic reconstruc- 
tions and palynological research (see OLALDE et al. 2002). 
This explains why many studies (e.g. TABERLET et al. 1998, 
PETIT et al. 2003) emphasized the role of glacial refugia in 
the Iberian, Apennine and Balkan peninsulas on the cur- 
rent distribution patterns of species. More recently, how- 
ever, phylogeographic analyses have accumulated evidence 
that many species may have found shelter in multiple 
small, isolated patches with milder climatic conditions (so- 
called ‘refugia within refugia, GÓMEZ & LUNT 2007). This 
scenario is more in line with the observed level of genetic 
variation and population sub-structuring in a wide variety 
of organisms, including amphibians, from the Apennine 
(e.g., CANESTRELLI & NASCETTI 2008), the Iberian (e.g., 
MARTINEZ-SOLANO et al. 2006) and the Balkan Peninsula 
(e.g., SOTIROPOULOS et al. 2007). Refugia are habitats with 
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space and time dimensions that operate on evolutionary 
time-scales and have facilitated the survival of species un- 
der changing conditions for millennia (KEPPEL et al. 2012). 
Consequently, identifying refugia, not only plays an im- 
portant role in understanding the evolutionary history of 
the world’s biota, but could even contribute to protecting it 
against future climate change (KEPPEL et al. 2012). Further- 
more, the high genetic diversity of populations that thrive 
in areas that acted as refugia during environmental insta- 
bilities, might represent the overall survival potential of a 
species, hence these populations may be the key for long- 
term species persistence (CONNER & HARTL 2004). 

In this study, we focus on the alpine salamander, Sala- 
mandra atra, a species endemic to the Alps, and the Di- 
narides (Fig. 1). Four subspecies have been described: 
(i) the completely melanistic (black) nominal subspecies, 
S. a. atra LAURENTI, 1768 from the Alps; two subspecies 
inhabiting the Italian Prealps, ie. (ii) S. a. aurorae TRE- 
VISIAN, 1982 and (iii) S. a. pasubiensis BONATO & STEIN- 
FARTZ, 2005, that are differentiated by the amount and pat- 
tern of yellow patches on the body (see BonaTo et al. 2018); 
and the fourth subspecies, also completely melanistic, (iv) 
S. a. prenjensis MIKŠIĆ, 1969, which consists of several iso- 
lated population fragments in the Dinarides. However, 
the validity of prenjensis has been debated for a long time 
(KLEWEN 1988, GROSSENBACHER 1994, GUEX & GROS- 
SENBACHER 2003, BONATO & STEINFARTZ 2005, DUBOIS 
& RAFFAELLI 2009, SPEYBROECK et al. 2010), although 
recent (evolutionary) studies support its separate taxo- 


nomic status (HELFER 2010, RAZPET et al. 2016, BONATO 
et al. 2018). Morphometric comparisons between selected 
populations of alpine salamanders suggest that prenjensis 
could have a larger head and jaws than the nominal sub- 
species (SuNjE et al. 2019), and, allegedly, it has a different 
arrangement of the palatal teeth (MIKŠIĆ 1969, KRIZMANIC 
1997). Although the contact zones between S. a. atra and 
S. a. prenjensis have not been defined, both subspecies are 
found in the Slovenian Alps (HELFER 2010, RAZPET et al. 
2016). From there, population fragments of alpine sala- 
manders are found across the Dinarides, more specifically 
at (from north to south): the mountain group of Trnovski 
gozd (Slovenia: Trnovski gozd, Hrušica, Menešija, Nanos) 
and Snežnik (Slovenia, RAZPET et al. 2016), Gorski Kotar, 
Žumberak, Istria and Plitvice (Croatia, JELIĆ et al. 2012), 
the mountains Prenj and Čvrsnica (Bosnia and Herzegovi- 
na, ŠUNJE & LELO 2010) and in the massive of Prokletije 
(Montenegro/Kosovo/Albania, KRIZMANIĆ 1997). Very re- 
cently, a new population fragment was discovered on Mt. 
Orjen (Montenegro, CIKOVAC & LJUBISAVLJEVIĆ 2020, 
Fig. 1). 

Dinaric samples of S. atra have been included in pre- 
vious phylogenetic studies on alpine salamanders: from 
the initial works of RIBERON et al. (2001, 2004) and Bo- 
NATO & STEINFARTZ (2005), up to the more recent studies 
of HELFER (2010) and BonaTo et al. (2018). However, the 
number of sampled specimens was invariably small and 
sampling was restricted to one or two locations in the Di- 
narides, hence the phylogenetic and phylogeographic rela- 
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Figure 1. Distribution of Salamandra atra with sampling sites (modified from Bonaro et al. 2018). Bordeaux polygons represent 
the known occurrence areas of S. atra, and the light blue lines are margins of the glaciers at their maximum extent during the last 
glacial period (Lgm, in the beggining of the Pleistocene, EHLERS & GIBBARD 2004); the grey color variation reflects altitude values, 
with higher altitudes being darker; the curve represents the border between the Dinarides and the Alps. A - Alpine, Prealpine, and 
North-Dinaric part of the distribution range. B - Distribution in the Central (Mt. Čvrsnica and Mt. Prenj) and Southern (Mt. Orjen 
and Prokletije Mts.) Dinarides. Full symbols are sampling sites from this study; empty symbols respresent sequences retrieved from 
public repositories (see text for details); code 22 (in grey) is both, sampled in this study and in Bonaro et al. 2018. 
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tionships among population fragments remained obscure. 
Hereby we present the first comprehensive study which 
includes most of the known population fragments of al- 
pine salamanders in the Dinarides. Using both microsat- 
ellite and mitochondrial (mt) DNA markers, we describe 
the evolutionary relationships, genetic structure and lev- 
els of differentiation among and within these isolated pop- 
ulation fragments. Further, as many taxa in the Balkans 
evolved following a ‘refugia within refugia scenario (e.g., 
Ichtyosaura alpestris, SOTIROPOULOS et al. 2007; Dinaro- 
mys bogdanovi, KRYŠTUFEK et al. 2007; Vipera ammodytes, 
URSENBACHER et al. 2008), and few experienced bottle- 
necks (e.g., Pinus nigra, NAYDENOV et al. 2017), we explore 
the possibility that such events occurred in the Dinaric 
population of S. atra. Considering the ‘refugia within refu- 
gia scenario, we expect a strong genetic subdivision among 
population fragments due to their isolation in separate re- 
fugia. This study contributes to a better understanding of 
the effects of past climatic oscillations on species confined 
to the mountainous habitats in the Dinarides, and provides 
relevant inputs to identify conservation units within the 
population of S. atra along the Dinaric arc. 


Material and methods 
Sampling 


Samples were collected in July and August 2016. We 
searched for individuals in suitable habitats focusing on 
the mountain areas in the Dinarides where S. atra has been 
confirmed: from north to south - Gorski Kotar (Samarske 
stijene and Vihora’ki put), Mt. Čvrsnica (Ploéno), Mt. 
Prenj (Podoti’, Kopilice and Zakantar) and Mt. Prokletije 
(Bogicevica and Gorazdevac). The sampling of individu- 
als in Kredarica (Julian Alps, sampling site: Aljažev dom - 
see Table 1) was done in 2010 and is included in this study 
for comparative purposes (unpublished data). Individuals 
were found on the soil or within shelters during favorable 
weather conditions (late night, early-morning or after in- 
tense rainfall). Small tail clips were taken from a total of 95 
individuals (Table 1). Clips were preserved in 96% ethanol. 


DNA extraction 


Total genomic DNA was extracted using the NucleoSpin® 
Tissue kit (Macherey-Nagel, Düren, Germany). DNA ex- 
traction from the samples of Kredarica was done in 2010 
using the protocol from TAGGART et al. (1992). All samples 
were brought to working concentrations by making 10-fold 
dilutions. 


mtDNA markers amplification and sequencing 
Two mitochondrial (mt) DNA markers were analyzed: the 


control region (D-loop) and cytochrome b gene (cob). Both 
regions were amplified and sequenced in two parts, using 


the Saat-H339-dloop/Saat-LPro-short (340 bp amplicon 
size) and L-Pro-ML/H-12S1-ML (820 bp) primers for D- 
loop, and the Sa_Cytb_12F/ Saat-H550 (560 bp) and Saat- 
L129/SaatCytb1046R (920 bp) primers for cob. All prim- 
ers were designed by Bonaro et al. (2018), except for L- 
Pro-ML/H-12S1-ML primer pair which was designed by 
STEINFARTZ et al. (2000). PCRs were performed in 25 ul 
volumes, following the protocols of BonaTo et al. (2018). 
The PCR products were sent to Macrogen Europe (Amster- 
dam, The Netherlands) for Sanger sequencing (using the 
forward PCR primers). 

The sequences were inspected in Jalview 2.9.0b2 (Wa- 
TERHOUSE et al. 2009) and manually corrected when 
needed; sequences were concatenated in Notepad (Micro- 
soft), aligned using ClustalW (THOMSON et al. 1994) and 
trimmed in Bioedit v5.09 (HALL 1999) to harmonize their 
length across the alignment. The protein-coding cob frag- 
ment was translated using the Translate tool on the ExPASy 
server (GASTEIGER et al. 2003) to check for open reading 
frames. The obtained sequences were submitted to Gen- 
Bank under the accession numbers (Acc. No.) MN255339 
— MN255350 and MN255326 - MN255337 (cob and D-loop 
respectively). Haplotype combinations (concatenated cob 
and D-loop) are given in Supplementary Table S1. 


Microsatellites genotyping 


Six nuclear microsatellite loci were amplified: SalE6, SalE7, 
SalE8, SalE12, SalE14 and SalE23 (STEINFARTZ et al. 2004) 
in two multiplex and one single reactions (Supplemen- 
tary Table S2). Each PCR reaction (total volume 10 ul) 
contained optimized concentrations of each primer pair 
(Supplementary Table S2), 1 x PCR buffer (20 mM Tris- 
HCl (pH 8.0), 100 mM KCl, 0.1 mM EDTA, 1 mM DTT, 
0.5% (v/v) Nonidet P4o, 0.5% (v/v) Tween 20, 50% (v/v) 
glycerol), 2 mM MgCl, 0.25 mM dNTPs, 1 U TaqNovaHs 
Polymerase (Blirt, DNA Gdansk, Poland) and 1 ul of DNA. 
All PCR reactions were done in Veriti Thermocycler (Ap- 
plied Biosystems) following the thermal protocols given in 
Supplementary Table S2. Genotyping was performed on 
the Genetic Analyzer 3500 Series (Applied Biosystems). 
PCR products from MIx 2 and the single reaction (SalE14) 
were pooled (1 ul each) and jointly genotyped. Allele scor- 
ing was done in GeneMapper v.5 using Liz 500 Size Stand- 
ard (Applied Biosystems). 


Phylogenetic and phylogeographic analysis 


DnaSP 6.0 (Rozas et al. 2017) was used to access the con- 
catenated mitochondrial sequence haplotypes (cob and D- 
loop). These haplotypes were analyzed together with 24 
other haplotypes (concatenated cob + D-loop) retrieved 
from published sequences (mainly Bonato et al. 2018, 
but also Vences et al. 2014, and STEINFARTZ et al. 2000 
for some cob sequences; see Supplementary Table S1) from 
118 individuals distributed over 19 sites across the species 
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Table 1. Dataset details and detected haplotypes (Hap) - Codes are as in Figure 1; codes from 17 to 27 are located in the Dinarides; 
bold codes are samples from this study, others are retrieved from public repositories (see text for details); code 22 is both: * - sampled 
in this study, and also ° - taken from Bonaro et al. (2018). The number of S. atra individuals included in the analyses of microsatellites 
(us) and concatenated - cob + D-loop - dataset (mt) is given in the respective columns; “/” - no data. Hap codes are as in Figure 2. At 
code 23, heteroplasmy (H6 + H12) was detected in one sample, hence four sequences were generated although only three individuals 
were sequenced. B&H - Bosnia and Herzegovina. 


Code Country/ Sampling area Sampling site us mt Hap (number of individuals) Lineage 

1 Italy / Breonie Pflerschtal / 6 H18 (1), H19 (1), H20 (4) S.a.atra 

2 Italy / Breonie Ridnauntal / 6 H19 (5), H20 (1) S.a.atra 

3 Italy / Dolomites Valle Fiscalina / 1 #H25(1) H25 

4 Italy / Orobie Val Tronella I 1  H29(1) Orobie 

5 Italy / Orobie Val Terzera / 8 H30 (4), H31 (4) Orobie 

6 Italy / Orobie Laghi Gemelli / 8  H27 (6), H28 (1), H32 (1) Orobie 

7 Italy / Orobie Pizzo della Presolana / 1 H27(1) Orobie 

8 Italy/ Orobie Val di Scalve / 9  H26 (1), H27 (8) Orobie 

9 Italy / Pasubio Val Fontana d'Oro / 13  #H14 (9), H15 (4) S.a.pasubiensis 
10 Italy / Sette Comuni Bosco del Dosso / 16 H21 (15), H23 (1) S.a.aurorae 
11 Italy / Sette Comuni Monte Fossetta / 1 H21 (1) S.a.aurorae 
12 Italy / Sette Comuni Val di Nos / 13  #H21 (11), H22 (2) S.a.aurorae 
13 Italy / Sette Comuni Val Postesina / i H21 (1) S.a.aurorae 
14 Italy / Belluno Schiara-Val dell’Ardo / 15 #H24 (15) H24 
15 Italy / Cansiglio Pian di Landro-Baldassare / 3 H24 (3) H24 
16 Slovenia / Kredarica Aljazev dom 8 8  H16 (1), H17 (7) S.a.atra 
17 Slovenia / Trnovski Gozd Trnovski Gozd / t 4H12(1) S.a.prenjensis 
18 Croatia / Gorski Kotar  Samarske stijene 7 4 H1 (1), H2 (1), H3 (1), H4 (1) S.a.prenjensis 
19 Croatia / Gorski Kotar Vihoraški put 10 5  H4 (4), H5 (1) S.a.prenjensis 
20 B&H / Čvrsnica Ploéno 20 9 H8 (8), H12 (1) S.a.prenjensis 
21 B&H/ Prenj Soplje / 7 H11 (2), H12 (4), H13 (1) S.a.prenjensis 
22 B&H/ Prenj Zakantar 7 3*+ 3° H9 (1°), H10 (1°) H12 (2*), H13 (1°+1*) S.a.prenjensis 
23 B&H / Prenj Kopilice 7 3 H6 (1), H12 (3) S.a.prenjensis 
24 B&H/ Prenj Podotis 11 4 H11 (1), H13 (3) S.a.prenjensis 
25 B&H/ Prenj Sedlo / 5 H11 (1), H12 (3), H13 (1) S.a.prenjensis 
26 Montenegro / Prokletije Bogicevica 22 9  H7(9) S.a.prenjensis 
27 Montenegro / Prokletije Gorazdevac 3 1 H7(1) S.a.prenjensis 


distribution range (Fig. 1, Table 1; the haplotypes of these 
individuals were provided by B. CRESTANELLO - Table S1). 
In total, 164 individuals were included in the phylogenetic 
and phylogeographic analysis (Table 1). 

To test for phylogenetic congruence, we performed 
a partition homogeneity test (FARRIS et al. 1995) in Paup 
4.0a (SWOFFORD 2002) with 1,000 replicates and heuristic 
search using the TBR branch swapping algorithm. Parti- 
tioning of the cob gene fragment (on 1°, 2" and 3 codon 
position) was done in DAMBE (X1a & XIE 2001), where- 
as the non-protein coding D-loop fragment was not par- 
titioned. DAMBE was also used to conduct a substitution 
saturation test (X1A et al. 2003) on the two mtDNA gene 
fragments. 

To build phylogenetic trees, sequences of D-loop and cob 
from Salamandra corsica Savi, 1838 (Acc. No.: MFo43388) 
and S. lanzai NASCETTI, ANDREONE, CAPULA & BULLINI, 
1988 (Acc. No.: MFo43391) were used as outgroups (follow- 
ing BonaTo et al. 2018). Phylogenetic relationships were 
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estimated using maximum likelihood (ML) and Bayesian 
inference (BI). ML analyses were conducted in RAxML v. 
8.2 (STAMATAKIS 2014) accessed through the CIPRES Sci- 
ence Getaway (MILLER et al. 2010) using the GTRGAM- 
MA model with 1,000 bootstraps. BI analyses were done 
in MrBayes v. 3.1.2 (RONQUIST & HUELSENBECK 2003) us- 
ing the best-fit nucleotide substitution model that was se- 
lected on the partitioned datasets and D-loop by the cor- 
rected Akaike Information Criterion (AICc) implemented 
in JModeltest 2.1.7 (DARRIBA et al. 2012) as suggested by 
BURNHAM & ANDERSON (2004). The analysis was run twice 
for 2 x 10° generations, with four parallel chains and sam- 
pled every 500 generations. 

The MCMC convergence was checked every 1,000 gen- 
erations by the output of standard deviation of split fre- 
quencies. The resulting tree was constructed from 6,000 
trees sampled from the posterior distribution, once the first 
2000 trees were excluded as burn-in (25%). Trees were vis- 
ualized and edited with FigTree v. 1.3.1 (RAMBAUT 2009). 
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Sequence divergence and geographical distribution of 
the concatenated haplotypes (invariable sites removed) 
were analyzed using unrooted median-joining networks 
(BANDELT et al. 1999) constructed with Poparr (http:// 
popart.otago.ac.nz). 


Genetic distances among lineages 


Arlequin 3.5 (EXCOFFIER & LISCHER 2010) was used to cal- 
culate the uncorrected pairwise genetic distances (pi; sim- 
ple p distances) for the concatenated mtDNA sequences 
among the lineages of S. atra inferred from the phyloge- 
netic analysis (see previously). The analysis was done un- 
der 10,000 permutations and a significance level of 0.05. 


Inferring the number of populations 


Using STRUCTURE 2.3.4 (PRITCHARD et al. 2010), the mi- 
crosatellites were analyzed to group the sampled individu- 
als of S. atra in clusters (K) which we define as populations. 
Following BonaTo et al. (2018), exploratory analyses, with K 
values between two and nine (the total number of sampling 
sites) were conducted with five simulations (100,000 iter- 
ations, 25,000 burn-in) implementing different combina- 
tions of ancestry models (admixture vs. no admixture) and 
allele frequency models (correlated vs. independent). The fi- 
nal analysis was performed under the admixture model and 
assuming correlated allele frequencies among populations 
(following RazpeT et al. 2016). The admixture model was 
chosen due to its flexibility (PRITCHARD et al. 2010), while 
the correlated allele frequency model was chosen because 
previous results showed genetic similarity between Dinar- 
ic population fragments (HELFER 2010, RAZPET et al. 2016, 
Bonaro et al. 2018). The analyses were run without using 
the information on sampling sites, and for K values from 2 
to 10, running 20 simulations (500,000 iterations, 125,000 
burn-in). The most likely number of clusters was estimat- 
ed by means of likelihood ratios (EVANNO et al. 2005) in 
STRUCTURE Harvester (EARL & VON HOLDT 2012). 

Using Genetix 4.05 (BELKHIR et al. 1996-2004), we per- 
formed a factorial correspondence analysis (FCA) on sam- 
pled individuals also to reveal potential population clusters 
and their genetic variation. 


Genetic variation within populations 


Microsatellite data were were checked for genotyping er- 
rors with Micro-checker 2.2.3 (VAN OOSTERHOUT et al. 
2004). Further analyses were performed in Arlequin. Devi- 
ation from Hardy-Weinberg equilibrium (HWE) was test- 
ed using the exact test (Guo & THOMPSON 1992). A link- 
age disequilibrium test was run between pairs of micro- 
satellite loci for each population with a likelihood-ratio test 
(EXCOFFIER & SLATKIN, 1998; with 10,000 permutations). 
Levels of significance (alpha = 0.01, see WIGGINTON et al. 


2005), were adjusted using the standard Bonferroni cor- 
rection. For each population we estimated the observed 
and expected heterozygosity, the mean number of alleles 
per locus and the inbreeding coefficient (Fis; tested using 
10,000 permutations of alleles between individuals with- 
in each population). The PopGenReport R package (ADA- 
MACK & GRUBER 2014) was used to assess the mean allelic 
richness of each population adjusted for sample size (us- 
ing a rarefaction procedure). The number of private alleles 
per population was obtained using the StrataG R package 
(ARCHER et al. 2016). 

mtDNA data: For the concatenated mtDNA sequenc- 
es, the number of (private) haplotypes, haplotype diversi- 
ty, number of polymorphic sites, and nucleotide diversity 
were computed in Arlequin. 


Genetic differentiation among populations 


The overall genetic variation was partitioned into within 
and among population components with an analysis of 
molecular variance (AMOVA, EXCOFFIER et al. 1992) us- 
ing Arlequin 3.5 (EXCOFFIER & LISCHER 2010); populations 
(as after STRUCTURE) were treated as one single group of 
data (one species, i.e. S. atra). The degree of among pop- 
ulation divergence was estimated by pairwise Fst (micro- 
satellite data, WEIR & COCKERHIAM 1984). Differentiation 
among mtDNA sequences was estimated by correcting for 
multiple hits with the method of TaMuRA & NEI (1993) be- 
cause this is the alternative for the HKY + I best fit evo- 
lutionary model (unavailable in Arlequin) that is detected 
within the cob gene (also following Bonaro et al. 2018). 
The significance of the differentiation values was assessed 
after 10,000 random permutations (EXCOFFIER et al. 1992) 
in Arlequin. 

A Mantel test was conducted using ade4 R package 
(Dray & DUFOUR 2007) with 9,999 permutations to evalu- 
ate whether the level of genetic differentiation among pop- 
ulations is correlated with their geographic distance. 


Demographic tests and bottleneck 


DnaSP 6.0 (Rozas et al. 2017) was used to compute 
Strobeck’s S, Fu and Li’s statistics and Tajima D test, for 
which we specified the coding (cob) and non-coding (D- 
loop) part of the sequence. The Tajima D test, combined 
with Fu and Li’s statistics, tests the hypothesis of range ex- 
pansions, while Strobeck’s statistics is relevant for the esti- 
mation of the number of existing haplotypes in relation to 
the number of recorded ones. 

The bottleneck hypothesis was tested for the inferred 
populations, using the one-tailed Wilcoxon test for hetero- 
zygosity excess implemented in Bottleneck 1.2.02 (Piry et 
al. 1999) under different mutation models (infinite allele - 
IAM, two-phased - TPM, stepwise mutation - SMM). The 
TPM model was tested with a variance of 14 and a propor- 
tion of SMM between 70 and 95% (Piry et al. 1999). 
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Results 
Phylogenetic and phylogeographic relationships 


To obtain a final consensus sequence including outgroups, 
we had to trim our concatenated sequences by 10 bp (the 
first two of cob and the last eight from the D-loop). The 
final sequence included 1,657 bp (cob = 962 bp, D-loop = 
695 bp). To achieve the same length, sequences of BONATO 
et al. (2018) had to be also trimmed, which caused sev- 
eral of their haplotypes to collapse with ours (Supplemen- 
tary Table S1). Finally, 32 distinct haplotypes were gener- 
ated (Supplementary Table S1). The partition homogene- 
ity test revealed no significant incongruence between cob 
and D-loop (p = 1). The substitution saturation test showed 
that cob partitions and D-loop are suitable for phylogenet- 
ic analysis. jModelTest revealed that cob evolution is best 
described under the models TPM2, HKY, HKY+I (1°, 2”¢ 
and 3" position respectively), while the D-loop evolution 
is under the TPM2UF + G model. The standard deviation 


S.lanzai 


0.63 /57 ------, H3 


0.99 / 68 ---- 
0.78 / 72 -- 
1.00/67 l. 


0.99 / 78 -- 


1.00/71 ...... ||]. -H18 3 
0.99 / 88 ~- |||. H19 
0.95 / 61 ------ H20 


0.90 / 58 
0.97/94 -----------F==4--- 


1.00 / 100 
0.98 / 64 


0.96 /62_ 
1.00/73 


1.00 / 99 H31 
Q 10 samples H32 
o 1sample 
0.004 


of the split frequencies after 2 x 10° generations was 0.009 
which is a good indication of convergence (RONQUIST et 
al. 2011). 

The clusters in the phylogenetic trees follow the ones in 
the haplotype network (Fig. 2) corresponding to those pre- 
viously revealed by BonaTo et al. (2018). The BI and ML 
revealed the same tree topology, suggesting a basal split 
between the Orobie population and all the others, and an 
unresolved polytomy between individuals from the Dolo- 
mites (H25), the populations in the Venetian Prealps (H24 
- Belluno and Cansiglio) and the clade uniting the pop- 
ulations assigned to the four subspecies (Fig. 2). The lat- 
ter clade separated populations from the Pasubio and Di- 
narides on one side, and the populations of Sette Comuni 
and the Alps on the other, although the monophyly of the 
first group is not well supported (branch values 0.84/49 - 
Fig. 2). In conclusion, the phylogenetic analysis revealed 
seven lineages, two of them being represented by single 
haplotypes (H24 and H25 respectively; Table 1, Fig. 2). 


Figure 2. Phylogenetic tree and Median-joining (MJ) network of the 32 haplotypes (H) based on the concatenated cob and D-loop 
markers in 164 individuals of Salamandra atra (46 samples from this study, and 118 from public repositories) from 27 sampling sites 
across the species range. Bayesian posterior probabilities (PP) and ML bootstrap values (in percentages from 1000 replications) are 
shown on the tree nodes separated by a slash (“/”); PP = 0.90 and ML values = 0.70 are bolded (considered high support values). In 
the MJ, the size of the haplotype (pie chart) is proportional to the number of samples in which it was found (as in the figure legend). 
The colors of haplotypes represent the sampling areas as in Table 1. The analysis revealed seven lineages : 1 - S. a. prenjensis (Dinaric 
clade, from north to south: black - Trnovski Gozd, yellow - Gorski Kotar, red - Čvrsnica, green - Prenj, pink - Prokletije); 2 - S. a. 
pasubiensis (Pasubio clade); 3 - S. a. atra (Alpine clade); 4 - S. a. aurorae (Sette Comuni clade); 5 - the Orobie clade potentially 
belonging to a new, yet undescribed subspecies; two lineages are represented by a single haplotype (H24 - brown, H25 - turquoise). 
The dashes along the branches of the MJ network represent the minimum number of substitutions between haplotypes that is also 
indicated by a (grey) number if different from one. Black circles in the MJ network represent missing haplotypes. 
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Table 2. Uncorrected p distances (pi) among Salamandra atra lineages; H24 - Venetian Prealps (Belluno and Cansiglio), H25 - Do- 
lomites, as in Table 1. Above diagonal: average pi between lineages; diagonal elements: average pi within lineages; below the diagonal 
are the p — values for estimates between the lineages (significant ones are in bold). 


Orobie S. atra prenjensis S. atra atra 
Orobie 1.28 27.29 23.05 
S. atra prenjensis 0.00 3.03 15.53 
S. atra atra 0.00 0.00 1.18 
S. atra aurorae 0.00 0.00 0.00 
S. atra pasubiensis 0.00 0.00 0.00 
H25 0.03 1.00 1.00 
H24 0.00 0.00 0.00 


S. atra aurorae S. atra pasubiensis H25 H24 
20.45 22.20 22.00 25.00 
16.19 14.53 16.22 16.44 

6.50 9.36 9.05 13.05 
0.90 11.76 11.45 15.45 
0.00 0.46 10.31 14.31 
0.06 0.08 0.00 12.00 
0.00 0.00 1.00 0.00 


Genetic distances among lineages 


According to the detected lineages from the previous 
analysis, we merged the sequences from the Dinarides 
(Trnovski Gozd, Gorski Kotar, Cvrsnica, Prenj and Prok- 
letije - sampling areas as in Table 1) in a single group 
named S. atra prenjensis; the sequences from Kredarica 
(Julian Alps, code 16 in Fig. 1) were merged with the se- 
quences from Breonie (codes 1 and 2 in Fig. 1) as these 
are part of the same lineage (S. a. atra, Table 1). Following 
the same principle we defined the other groups represent- 
ing the inferred lineages (S. a. aurorae, S. a. pasubiensis, 
Venetian Prealps - H24, Dolomites - H25, Orobie, Ta- 
ble 1). The analysis of p-distances showed that all lineages 
are significantly distant among each other (Table 2) ex- 
cept H25 which differs significantly only from the Oro- 
bian lineage (p = 0.03, Table 2). However, taking into con- 
sideration that only one individual represents lineage H25 
in the Dolomites (code 3, Table 1) the estimated distance 
between this and other lineages (Table 2) must be taken 
with precaution. 


Summary statistics for samples originating 
from this study 


A total of 95 individuals of S. atra were typed for six mi- 
crosatellite loci (Table 1). All loci were polymorphic, with 
seven to 12 alleles per locus (mean = 8.7, Supplementary 
Table S2). No null alleles, allelic dropout, or stuttering were 
detected for any locus. Raw microsatellite data are given in 
Supplementary Table S3. 

The alignment of the concatenated (consensus) mtDNA 
sequences had a length of 1667 bp (964 bp cob + 703 bp 
D-loop) and comprised 46 individuals from nine sampling 
sites (Table 1). Within the D-loop, heteroplasmy was ob- 
served (506 bp, C-T) in one individual from Prenj (sam- 
pled in Kopilice), therefore the two alternative copies of 
this sequence were kept (Accession numbers: MN255332 
and MN255334). The Dinaric samples involved 13 haplo- 
types (Fig. 2), with 28 polymorphic nucleotide positions 
(1.6%; including nine transitions, four transversions, one 
indel and 14 substitutions [0.8%]). 


Inferring the number of populations 


All exploratory model combinations performed in STRUC- 
TURE proposed five clusters (Fig. 3, Supplementary Fig. S1). 
Although each cluster (population) mainly corresponds to a 
separate mountain area (from north to south: Kredarica (Ju- 
lian Alps), Gorski Kotar, Mt. Čvrsnica, Mt. Prenj, Mt. Prok- 
letije), the clusters are not completely differentiated because 
some individuals in each population, resemble individu- 
als from others with a high probability (especially noted in 
Gorski Kotar, Čvrsnica [1 sample] and Prokletije, Fig. 3). 

The first FCA axis revealed two main groups: the first 
is more variable and contains only the samples from the 
Julian Alps (Kredarica), and the second groups all the Di- 
naric samples. Within the second group, four fairly distinct 
clusters can be recognized, each corresponding to a sepa- 
rate mountain area (Fig. 4). The combination of the three 
FCA axes accounted for 24.02% of the total variation. No 
clear substructure was detected among individuals of the 
same mountain area (clustering of individuals from differ- 
ent sampling sites within the same mountain area, Supple- 
mentary Fig. S2). 


Genetic variation within populations 


Estimated diversity indices for both microsatellites and 
mtDNA are shown in Table 3. Only the population of Prok- 
letije deviated from HWE for the locus SalE7 (p = 0.0001, 
after Bonferroni correction). No deviations from linkage 
equilibrium were detected. 

The observed microsatellite heterozygosity for Dinar- 
ic specimens varied between 0.31 and 0.39, but was nearly 
twice as high in the Alpine population (Kredarica; Table 3). 
The highest and only significant inbreeding coefficient (Fis) 
was recorded in Gorski Kotar (0.19, p = 0.009, Table 3). 

The mtDNA diversity (pS, Hap, pHap, Hd in Table 3) 
was higher in Prenj and Gorski Kotar, and lower in Prok- 
letije, Cvrsnica and Kredarica. Only one (private) haplo- 
type was found in the population of Prokletije. The high- 
est number of haplotypes (N = 5, all private) was observed 
in Gorski Kotar, but the highest nucleotide diversity - in 
Prenj (Nd = 0.18, Table 3). 
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Table 3. Genetic diversity from microsatellite data (us - in grey) and mtDNA haplotypes (concatenated cob and D-loop - in white) in 
the populations (Pop) of Salamandra atra from this study; N - number of individuals (total microsatellites - N = 95, total mtDNA - 
N = 46), Ho - observed heterozygosity, He - expected heterozygosity, mA - mean number of alleles, R - allelic richness, pA - number 
of private alleles, Fis - inbreeding coefficient (bold if significant); pS - number of polymorphic sites, Hap - number of haplotypes, 
pHap - number of private haplotypes, Hd - haplotype diversity, Nd - nucleotide diversity (averaged over the two gene fragments). 
Numbers in brackets are standard deviations. 


Pop N Ho He mA R pA Fis N ps Hap pHap Hd Nd (%) 
Kredarica 8 0.58 (0.25) 0.58 (0.19) 3.83 (0.98) 3.60 10 -0.002 8 1 2 2 0.25 (0.18) 0.02 (0.02) 


Gorski Kotar 17 0.31 (0.21) 0.39 (0.26) 3.83 (1.47) 2.93 7 0.19 9 5 5 5 0.72 (0.16) 0.05 (0.05) 


Prenj 
ren) 25 0.39 (0.28) 0.36 (0.28) SR 2.78 4 -0.07 10 7 4 3 0.71 (0.10) 0.18 (0.11) 
Cvrsnica 
ZOMOS IN (OFS2) MO SON(O27) o (e754 Ones -0.03 9 2 2 1 0.22 (0.17) 0.03 (0.03) 
Prokletije 25 OS(s) O27 (O23) aes) B72 5 0.09 10 0 1 1 0 0 


Figure 3. Results of the STRUCTURE assignment of 95 individuals of Salamandra atra (from this study) based on six microsatellite 
loci; clusters (populations) are divided by tick black lines as follows: 1 - Kredarica, 2 - Gorski Kotar (. - Vihoraski put, .. - Samarske 
stijene), 3 - Čvrsnica, 4 — Prenj (. - Kopilice, .. - Podotiš, ... - Zakantar), 5 - Prokletije (no dot - Bogi¢evica, . - Gorazdevac,). Com- 
pared to the Kredarica population (population 1 - in the Julian Alps), the Dinaric ones (populations - 2-5) are less differentiated, as 
some individuals in each population resemble individuals from others. 


Š K T 
N p20 gogo 
i a 
Si ° 2 a B E gak 
N a oF Pon a0 of a a 
3 CHE a 
-1 erry 
. " Axis 1 (9.78%) ° 5 j ° Axis 2 (7.54%) i 


Figure 4. Factorial Correspondence Analysis (FCA) ordination along the first three axes of 95 individuals of Salamandra atra from four 
sampling areas in the Dinarides (following Table 1) and one in the Julian Alps (Kredarica), based on six microsatellite loci. Each label 
corresponds to a sampling area (from north to south): blue - Kredarica, yellow - Gorski Kotar, red - Čvrsnica, green - Prenj, pink - 
Prokletije. Axis 1 separates the Dinaric individuals on one side, and the ones from Kredarica on the other; Axis 3 separates Gorski 
Kotar (Northern Dinarides) from the other sampling areas in the Dinarides (Central and Southern - see also Supplementary Fig. S2). 
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Genetic differentiation among populations 


With the microsatellite data, the among-population dif- 
ferentiation accounted for a significant 35.09% of the total 
variation (Amova, p < 0.001). Pairwise Fst values ranged 
from 0.28 (between Gorski Kotar and Prenj) to 0.46 (be- 
tween Prokletije and Cvrsnica) and were all significant 
(p < 0.001, Table 4). 

With the mtDNA data, the among-population differen- 
tiation accounted for a significant 87% of the total varia- 
tion (AMova, p < 0.001). Pairwise mtDNA values (@ val- 
ues - TAMURA & NEI 1993) between populations ranged 
from 0.43 to 0.99, and were all significant (p < 0.001, Ta- 
ble 4). The strongest differentiation was observed between 
Kredarica and all other populations; within the Dinar- 
ides, mtDNA showed the highest differentiation between 
Čvrsnica and Prokletije (0.94) and the lowest between 
Prenj and Gorski Kotar (0.43; Table 4) which is concordant 
to previously reported pairwise Fst values. 

The Mantel tests revealed no correlation between genet- 
ic differentiations and geographic distances (p 
0.95 Panna = 0-6). 


microsatellites 


Demographic tests and bottleneck 


The Tajima D test could not reject neutral mtDNA se- 
quence evolution (coding region: D = -0.65, p > 0.10; syn- 
onymous sites: D = -0.16, p > 0.10; nonsynonymous sites: 
D = -0.24, p > 0.10; silent sites: D = -0.24, p > 0.10), and 
thus did not provide evidence of range expansions. The 
results of the Tajima D test were in line with those of Fu 
and Li's statistic, which neither could reject neutrality (D = 
-1.27, p > 0.10; F = -1.24, p > 0.10). Strobeck’s S statistic pre- 
dicted a higher number of haplotypes than it has been ob- 
served (S =1, p = 0.00). 

The bottleneck tests were unable to demonstrate that 
such events occurred in our populations (IAM: p = 0.50, 
Smm: p 2 0.96, Tpm: p 2 0.44). 


Discussion 
Phylogenetic and phylogeographic relationships - 
taxonomic implications 


All studied populations in the Dinarides belong to a 
unique evolutionary lineage (Fig. 2). Individuals from the 
Julian Alps (Kredarica) belong to S. a. atra (Table 1, Fig. 2) 
as it could have been anticipated from the high differen- 
tiation values between this and all the Dinaric populations 
(mtDNA, Table 4), and its high heterozygosity (Table 3) 
that are characteristic to this lineage when compared to 
all the others (HELFER 2010, RAZPET et al. 2016, BONATO 
et al. 2018). This study confirmed already observed trends 
considering the evolutionary separation of the Dinaric 
population and the delineation of the commonly accept- 
ed subspecies: S. a. aurorae, S. a. pasubiensis and S. a. atra 
(HELFER 2010, BonaTo et al, 2018). Moreover, the study 


Table 4. Pairwise divergence between populations of S. atra from 
this study; pairwise Fst for microsatellite loci is above the di- 
agonal, and differentiation among mtDNA loci (cob + D-loop) 
is below the diagonal (ọ values - TAMURA & NEI 1993, see text 
for details); all values (for both microsatellite and mtDNA) were 
significant (p < 0.001). 


Kredarica eo Prenj Cvrsnica Prokletije 
otar 
Kredarica / 0.30 0.36 0.33 0.33 
Gorski Kotar 0.96 / 0.28 0.41 0.30 
Prenj 0.88 0.43 / 0.35 0.35 
Čvrsnica 0.98 0.78 0.53 / 0.46 
Prokletije 0.99 0.88 0.61 0.94 / 


also confirmed the evolutionary separation of the Orobi- 
an population of alpine salamanders and the populations 
from the Venetian Prealps and Dolomites (H24 and H25 
respectively, Table 1, see also HELFER 2010, HELFER et al. 
2011 and BonaTo et al. 2018). However, as our phylogenetic 
results were inferred from mtDNA data, conclusions need 
further corroboration, since in S. atra, nuclear markers 
may generate different tree topologies (see BONATO et al. 
2018, RIBERON et al. 2004 vs. HELFER 2010). By reconfirm- 
ing the separate evolutionary history of our Dinaric speci- 
mens and by presenting their significant genetic distance 
from other lineages (Table 2), we justify the subspecific 
taxon S. a. prenjensis MIKŠIĆ, 1969 as already proposed by 
previous authors (HELFER 2010, RAZPET et al. 2016, BONA- 
TO et al. 2018, ŠUNJE et al. 2019). 

In order to fill the remaining gaps in our knowledge 
about the complex evolutionary history of S. atra, we high- 
light the importance to study other population fragments 
that are not included in this study; primarily the ones in 
the Prealps and Dinarides (Fig. 1, see also HELFER 2010). 


Inference of populations 


Both STRUCTURE and FCA showed that the specimens 
from this study belong to five populations that mainly cor- 
respond to the sampled mountain areas. Populations are 
well differentiated but not completely, which is especially 
visible within the Dinaric clusters where in nearly all popu- 
lations some individuals are genetically more similar to in- 
dividuals from other populations than to members of their 
own population (Fig. 3). Considering the large geographic 
distance among the population fragments and the strong 
geographic barriers between them (Fig. 1), it is improba- 
ble that the incomplete population differentiation is due to 
current gene flow. Indeed, a lack of gene flow is confirmed 
by the high and significant differentiation values (Ta- 
ble 4 and Amova). Instead, the incomplete differentiation 
among the Dinaric specimens may be due to a past con- 
nection between the populations and/or a common origin. 

The population of Kredarica was the most homogenous 
(Fig. 3) and most distantly related (Fig. 4) compared to all 
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Dinaric specimens (S. a. prenjensis lineage - Table 1, Fig. 2) 
which is explained by the fact that it belongs to another lin- 
eage (S. a. atra — Table 1, Fig. 2). 


Distribution of genetic diversity 


The observed heterozygosity of the Prokletije population 
(Ho = 0.33) is twice as high as the Ho reported by Raz- 
PET et al. (2016; Ho = 0.17) probably because the latter 
study involved a smaller number of individuals from this 
location, two of which were intrauterine larvae related to 
a third sample (mother). The observed heterozygosity in 
Prenj was in line with previous observations (RAzPET et al. 
2016, BONATO et al. 2018: Ho = 0.36 and 0.37 vs. Ho = 0.39 
in Table 3). The microsatellite genetic diversity indices were 
similar in all the Dinaric populations, whereas for mtDNA, 
the levels of genetic diversity are higher in Gorski Kotar 
and Prenj, and lower in Cvrsnica and Prokletije (Table 2). 
On the other hand, the heterozygosity indices for Kredari- 
ca, were almost twice as high as in the Dinaric populations 
(despite a considerably smaller sample size) but, the mt- 
DNA diversity was lower than in Cvrsnica and Prokletije 
(probably due to a considerable smaller sample size - Ta- 
bles 1 and 3). A high heterozygosity in Kredarica (both ob- 
served and expected) suggests gene flow that may be a con- 
sequence of a “leading edge” range expansion characteris- 
tic for other populations of S. atra atra across the north- 
ern and central Alps (HELFER 2010). Just like RAzPET et al. 
(2016) we were unable to detect any signs of bottlenecks. 


‘Refugia within refugia’: 
a plausible historical scenario for the prenjensis lineage 


The ‘refugia within refugia hypothesis describes survival 
patterns in multiple glacial refugia, each within a larger 
refugial area (GOMEZ & LUNT 2007). The Balkan refugi- 
al area was mostly free of ice during the last glacial maxi- 
ma (LGM, Fig. 1) but many glaciers persisted during oth- 
er phases of the Pleistocene at higher altitudes (generally 
above 1,300 m a.s.l., CVIJIĆ 1899, HUGHES et al. 2011, LEPI- 
RICA 2013, ZEBRE & STEPISNIK 2015), which explains why 
multiple refugia must have characterized the Balkans. Our 
results suggest that during the Pleistocene, Dinaric popu- 
lations of alpine salamanders thrived in at least two refu- 
gia: an ancestral in Mt. Prenj and, a more recent, in Gor- 
ski Kotar. The diversity pattern in Prenj (the highest Nd, 
high Hd with many pHap, Table 3) and its central position 
in the MJ network (Fig. 2), suggest that it represented the 
main diversification center of alpine salamanders (follow- 
ing AVISE 2000). The fact that the lowest mtDNA differen- 
tiation is found between Prenj and all other populations 
(Table 4) supports this idea. Moreover, RAzPET et al. (2016) 
showed that Mt. Prenj is home to the oldest populations of 
alpine salamanders in the Dinarides (according the esti- 
mated time since decline). Gorski Kotar was another di- 
versification center of alpine salamanders in the Dinarides, 
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as is suggested by its high mtDNA diversity, which resem- 
blances the pattern found at Prenj (Table 3). However, the 
diversification at Gorski Kotar appears to have a more re- 
cent origin given its low nucleotide diversity and its star- 
like network of very similar haplotypes (Fig. 2, following 
AVISE 2000). Hence, Prenj and Gorski Kotar, interpreted 
as alpine salamander refugia, might not have co-existed. 
Furthermore, the central haplotype in Prenj (H12) is the 
most frequent in the ‘prenjensis’ clade (Fig. 2, Table 1), sug- 
gesting that it is a very old haplotype that was probably dis- 
tributed over the entire Dinaric range and persisted in big 
populations but was lost in others during fragmentations. 
This ancestral haplotype is separated by only one single 
mutational step from the central haplotype in Gorski Kotar 
(H4, Fig. 2), suggesting that alpine salamanders at Gorski 
Kotar originate from Prenj. Nevertheless, the many private 
haplotypes found in Prenj and Gorski Kotar are also evi- 
dence of a long-term population establishment (see Hay- 
DAR 2011), which supports our conclusion that these areas 
acted as refugia for S. atra. 

In contrast, the lower mtDNA diversity and fewer pri- 
vate haplotypes in Cvrsnica and Prokletije (Table 3), sug- 
gest that these populations were established more recently 
(following AUSTERLITZ et al. 1997 and ExCOFFIER & PETIT 
2009). Considering Čvrsnica, these results are rather un- 
expected, as this mountain is neighboring Prenj. Although 
several areas of both mountains were glaciated in the mid- 
dle Pleistocene (after the LGM; Mt. Čvrsnica - above 
1,200 m a.s.l. [STEPIŠNIK et al. 2016], Mt. Prenj - above 
1500 m [LEPIRICA 2008]) the population of Prenj has a 
considerably higher mtDNA diversity than Cvrsnica. We 
suspect that the relatively larger glacier on Cvrsnica (ca. 
200 km’, see also MILIČEVIĆ 2013) hindered the popula- 
tion stability more than did the glacier on Mt. Prenj. This 
is especially plausible given that the survival of individuals 
on Mt. Čvrsnica was thwarted by several barriers: on the 
east, by the deep canyon of the river Neretva (at the foothill 
of Mt. Čvrsnica), on the west — by the glaciers Svinjača and 
Dugo polje (also on Mt. Čvrsnica, STEPIŠNIK et al. 2016), 
and on the south by the glacier of Mt. Cabulja (south bor- 
der of Mt. Cvrsnica, see PRSKALO 2008). Moreover, the 
highest differentiation (both Fst and mtDNA) found be- 
tween Čvrsnica and Prokletije (the southernmost popula- 
tion; Table 4) confirms that there is no dispersal from Mt. 
Čvrsnica southwards. The missing haplotype link between 
Čvrsnica and both, Gorski Kotar and Prenj (separated by a 
single mutational step, Fig. 4) suggests that the stabilization 
of Cvrsnica occurred from another population (not sam- 
pled/not existing) from which dispersal was probably bi- 
directional, i.e from north (Gorski Kotar) to south (Prenj) 
and vice versa (Fig. 2). Although speculative, the HWE de- 
viation at SalE7 may be a hint that this area was colonized, 
as founding events are likely to create persistent non-equi- 
librium structures (BOILEAU & HEBERT 1991). Individuals 
that were part of this founding event likely originated from 
Prenj since the only haplotype registered in Prokletije (H7, 
Fig. 4) is separated from the Prenj haplotype by only one 
mutational step (H6, Fig. 2 and Table 1). We conclude that 
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our results are in line with the ‘refugia within refugia hy- 
pothesis that was suggested for several other species in the 
Balkans (KRYŠTUFEK et al. 2007, SOTIROPOULOS et al. 2007, 
URSENBACHER et al. 2008, SuRINA et al. 2011). 


Conservation insights for the prenjensis lineage 


A major conservation task is to define the proper distri- 
bution range of prenjensis. Strobeck’s S statistic and the 
newly discovered population on Mt. Orjen (CIKOVAC & 
LJUBISAVLJEVIC 2020) suggest that there are populations 
of alpine salamanders in the Dinarides that are yet to be 
discovered (see also KLEWEN 1988). Therefore, field inves- 
tigations of potential occurrence areas (see in CIKOVAC & 
LJUBISAVLJEVIC 2020) are needed to clarify the distribu- 
tion range of prenjensis. The northernmost locations where 
prenjensis is reported are in the Southern Prealps (Koba- 
rid - Krn, HELFER 2010 and Loibl Pass, Sunyz et al. 2019). 
Since these two areas are relatively close to our sampling 
site of Kredarica (ca. 20 and 30 km of air distance respec- 
tively), several contact zones between S. a. atra and S. a. 
prenjensis may be present in the wider areas connecting the 
three sites. The identification of contact zones is, not only 
important from an evolutionary perspective and the study 
of speciation processes (WOLLENBERG VALERO et al. 2019), 
but it is relevant in a conservation context (ALLENDORE et 
al. 2001), since conservation efforts should be preferably 
oriented towards populations where there is no risk of ge- 
netic swamping with other lineages. Genetic swamping 
may cause outbreeding depression (FRANKHAM et al. 2011), 
that, at it’s extreme form, could be detrimental for the pop- 
ulation (RHYMER & SIMBERLOFF 1996). The levels of genet- 
ic diversity of the Dinaric populations in this study suggest 
that these are not of immediate conservation concern, al- 
though the significant level of inbreeding in the population 
of Gorski Kotar (Fis, Table 3), signals a potential threat to 
its health and must be monitored. The significant differen- 
tiation among populations, indicates that each should be 
treated as a separate conservation unit. As the population 
of Prenj harbors substantial genetic diversity and no signs 
of inbreeding (Table 3), we believe that it is in conservation 
interest to focus the efforts on its preservation on the long 
term. Another reason to focus conservation efforts on the 
Prenj population is that Prenj was a refugium for S. atra 
during past environmental perturbances, so that it might 
also offer the best chances for survival under upcoming cli- 
matic change (for the rationale see KEPPEL et al. 2012). 
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Supplementary data 


The following data are available online: 


Supplementary Figure S1. STRUCTURE exploratory analysis 
using different combinations of ancestry and allele frequency 
models 


Supplementary Figure S2. Factorial Correspondence Analysis 
(FCA) ordination along the first three axes of 67 individuals of 
Salamandra atra from the Dinarides based on six microsatellite 
loci. 

Supplementary Table S1. Accession numbers (Acc. No) of the 
sequences (cob and D-loop) from public repositories (GenBank) 
used in this study. 


Supplementary Table S2. PCR conditions and characteristics of 
the microsatellite loci. 


